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Introduction for "An Introduction to Wavelet Analysis" 


In this chapter, we give an overview on multiresolution analysis, wavelet 
series and wavelet estimators in the classical setting. By ‘classical’ or ‘first- 
generation’ wavelets, we mean wavelets that were constructed initially to 
analyze signals observed at equispaced design points and having a sample 
size which is a power of two. The ‘second-generation' wavelet basis 
presented in the subsequent chapters will release these two constraints. 


If one wants to analyze a function of time with a series expansion, the first 
idea that comes probably into one's mind is to use a Fourier series, i.e. 
decompose the function into sine and cosine at different frequencies. In this 
process, we hope that only a few coefficients in the series will carry most of 
the information about the signal. Certain smooth functions have such an 
“economical' Fourier expansion. However, for most functions, a good 
Fourier series approximation requires numerous sine and cosine basis 
functions. Indeed, the sine functions have a precise frequency but are not 
localized in time, hence a localized information in the signal like a 
discontinuity will affect all the coefficients of the series. This drawback 
lead the researchers to look for more efficient bases, that is, bases which are 
localized both in time and in frequency. We will see in Multiresolution 
analysis and wavelets that a wavelet basis offers this property. 


This chapter is structured as follows. We begin by recalling some notations 
in Function spaces: notion and notations. Next Multiresolution analysis and 
wavelets introduces the multiresolution analysis, the wavelet functions, and 
gives some simple examples of wavelet bases. Fast wavelet transform 
explains how to decompose a signal using the wavelet transform. Such 
wavelet transforms, also called ‘decimated’, lack the property of translation- 
invariance. Non-decimated wavelet transform presents a widely used trick 
to make a wavelet transform translation-invariant. Since the main goal of a 
wavelet series is to provide a good approximation of a function belonging 
to a given space, Approximation of Functions introduces some fundamental 
notions to measure the quality of such approximation. Finally, 
Nonparametric regression with wavelets presents how to construct a 
nonparametric regression estimator using wavelets. First, the classical case 
of equally spaced design is considered. In the last part of Nonparametric 


regression with wavelets, we review some existing methods that deal with 
irregular designs. 


Function spaces: notion and notations 


A Hilbert space is a complete normed space whose norm is indexed by an 
inner (or scalar) product. 


Two disjoint subspaces A and B of a space S form a direct sum 
decomposition of S if every element of S can be written uniquely as a sum 
of an element of A and an element of B. The notation S = A @ B is then 
used. 


A measurable function f belongs to the Lebesgue space 
L,(R), 1 <p< wif 


Equation: 
+00 1/p 
itt, = ( f e))") <0. 


An example of a Hilbert space is the Lebesgue space D2 (IR) of measurable 
and square integrable functions. Indeed, the norm ||-||, is induced by the 
scalar product 

Equation: 


Fg / f (a)g(a)de , 


where g(x) denotes the complex conjugate of g(x). Two functions are said 
to be orthogonal in D2 (R) if their inner product is zero. 


The Lebesgue measure can be replaced by a more general measure jJ, 
leading to the weighted space Le (1), which has as inner product 
Equation: 


(faa, = / f (e)g(a)dp (e) 


and which contains the functions that have a finite norm 


IF llan = VF FY au Eto en 


A countable subset { f;,,} of functions belonging to a Hilbert space is a Riesz 
basis if every element f of the space can be written uniquely as 

f = >, crf, and if positive constants A and B exist such that 

Equation: 


Allfllz < d= leel? < BI FI; - 
k 


A Riesz basis is an orthogonal basis if the f;, are mutually orthogonal. In 
this case, A= B= 1. 


Multiresolution analysis and wavelets 


Definition of subspaces and of scaling functions 


A natural way to introduce wavelets is through the multiresolution analysis. Given a function f € D2 (R),a 
multiresolution of L» (IR) will provide us with a sequence of spaces V;, Vj+1, ... such that the projections of f onto 
these spaces give finer and finer approximations (as 7 —> oo) of the function f. 


Note: (Multiresolution analysis (MRA) in the first generation) A multiresolution analysis of L2 (IR) is defined 
as a sequence of closed subspaces V; C D2 (R), 7 € Z with the following properties: 


1. Equation: 


= CVYACVCVAC... 


2. The spaces V; satisfy 
Equation: 


U V; is dense in L.(R) and ‘al V; = {0} . 


jEZ jeZ 


3.1f f(z) Ee Vo, f (2/2) € Vj, i.e. the spaces V; are scaled versions of the central space Vo. 
4. If f € Vo, f (.—k) € Vo,k € Z, that is, Vo (and hence all the V;) is invariant under translation. 
5. There exists ¢ € Vo such that {¢(a — k);k € Z} is a Riesz basis in Vo. 


We will call ‘level’ of a MRA one of the subspaces V;. From [link], it follows that, for fixed j, the set 

{ 5h (x) = 29/24 (2?¢ = k) ;ke Zz} of scaled and translated versions of ¢ is a Riesz basis for V;. Since 
@ € Vo C Vj, we can express ¢ as a linear combination of {1}: 

Equation: 


(a) =) hadie(w) = V2.7 had (2a — k) - 


keZ keZ 


[link] is called the two-scale equation or refinement equation. It is a fundamental equation in MRA since it tells 
us how to go from a fine levelV; to a coarser levelVo. The function ¢ is called the scaling function. 


As said before, the spaces V; will be used to approximate general functions. This will be done by defining 
appropriate projections onto these spaces. Since the union of all the V; is dense in D2 (IR), we are guaranteed that 
any given function of D» (IR) can be approximated arbitrarily close by such projections. As an example, define the 
space V; as 

Equation: 


Vj = {f € Lo (R); Vk € Z, flask 2v(rery = constant} 


Then the scaling function ¢ (a) = 1(9,1) (x), called the Haar scaling function, generates by translation and 
dilatation a MRA for the sequence of spaces WA j € Z} defined in [link], see [link], [link]. 


The detail space and the wavelet function 


Rather than considering all the nested spaces V;, it would be more efficient to code only the information needed to 
go from V; to V;,1. Hence we consider the space W; which complements V; in Vj1 : 
Equation: 


Vi = Vi @ W; . 
The space W; is not necessarily orthogonal to V;, but it always contains the detail information needed to go from 
an approximation at resolution 7 to an approximation at resolution 7 + 1. Consequently, by using recursively the 


equation [link], we have for any jo € Z, the decomposition 
Equation: 


Lz (R) = Vj, © ©%2,,Wj - 


J=Jo 


With the notational convention that W;,_; := Vj,, we call the sequence 
{W3} 555-1 a multiscale decomposition (MSD). 
We call w a wavelet function whenever the set {W(x — k); k € Z} is a Riesz basis of Wo. Since Wo C Vi, there 


also exist a refinement equation for w, similarly to [link]: 
Equation: 


b(2) = V2S— ond (2a — k) . 


The collection of wavelet functions {Wik = 25/24) (27x — k) *kKEZ,IE Zz} is then a Riesz basis for Lz (R). One 
of the main features of the wavelet functions is that they possess a certain number of vanishing moments. 


Note:A wavelet function 7)(x) has Nvanishing moments if [ ~(x)x?dx = 0, p=0,...,N —1. 


We now mention two interesting cases of wavelet bases. 


Orthogonal bases 


In an orthogonal multiresolution analysis, the spaces W;; are defined as the orthogonal complement of V; in 
Vj41. The following theorem tells us one of the main advantages of such a MRA. 


Note:([link], Theorem 5.1.1) If a sequence of closed subspaces (V;), «z mL, (R) satisfies [link], and if, in 


addition, {¢(2 — k), k € Z} is an orthogonal basis for Vo, then there exists one function ~(a) such that 
{y(a — k); k © Z} forms an orthogonal basis for the orthogonal complement Wg of Vo in V4. 


An immediate consequence of [link] is that {wj,,k € Z} constitutes an orthogonal basis for the orthogonal 
complement W;, of V; in V;,1. In this section, let Y; (resp. 2;) be the orthogonal projection operator onto V; 
(resp. W;). The orthogonal expansion 

Equation: 


f = Pyit+)> 2f 


I=Jo 
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tells us that a first, coarse approximation of f in V;, is further refined with the projection of f onto the detail 
spaces W;. 


[link] shows two examples of orthogonal wavelet functions. The first is the Haar wavelet, associated to the Haar 
scaling function defined in !'Definition of subspaces V_j_and of scaling functions". 
Equation: 


Pla) = 2-0? (p49 (Qe — 1) — pH (x) = 1 p23) (@) — Lory (2) - 


The Haar wavelet has only one vanishing moment and consequently is optimal only to represent functions having a 
low degree of regularity, like, for example, 8—Hd6lder functions with 0 < 6 < 1. 


Daubechies constructed in [link], [link] compactly supported wavelets which have more than one vanishing 
moment. Compactly supported wavelets are desirable from a numerical point of view, while having more than one 
vanishing moment allows to reconstruct exactly polynomials of higher order. These wavelets cannot, in general, be 
written in a closed analytic form. However, their graph can be computed with arbitrarily high precision using a 
subdivision scheme algorithm. [link](b) represents the Daubechies Least Asymmetric wavelet with N = 4 
vanishing moments. 


Some orthogonal basis functions: (a) the Haar wavelet function bases with N = 1 vanishing moments, (b) the 
Least Asymmetric wavelet function of Daubechies [link], [link], with N = 4 vanishing moments. 
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This figure also illustrates the reason behind the name 
wave < tl: sin cewave < tsarefunctionswithacerta € vmberofvanish € gmoments, theyhavetheshapeo,. 
little wave' or ‘wavelet’. 


Biorthogonal bases 


Having an orthogonal MRA puts strong constraints on the construction of a wavelet basis. For example, the Haar 
wavelet is the only real-valued function which is compactly supported and symmetric. However, if we relax 
orthogonality for biorthogonality, then it becomes possible to have real-valued wavelet bases of fixed but 
arbitrary high order (see Definition 1 from Approximation of Functions) which are symmetric and compactly 


supported [link]. In a biorthogonal setting, a dual scaling function b and a dual wavelet function wy exist. They 
generate a dual MRA with subspaces V; and complement spaces W; such that 
Equation: 


V; | W; and V; 1 W;. 


In other words, 
Equation: 


(G,u(-—k))=0 and (6,8(-—h)) =0 


Moreover, the dual functions also have to satisfy 
Equation: 
(5,6 = k)) = 6k,0 and (0, Hl = k)) = 0k0 » 


where 4,9 is the Kronecker symbol. By construction, the dual scaling and wavelet functions satisfy a refinement 
equation, similarly to the equations [link] and [link]. 


In this work, we use the following convention: the dual MSD will be used to decompose a function (or a signal), 
while the original, or primal MSD reconstructs the function. This yields the following representation of a function 
f € L2(R) 

Equation: 


f («) = LaF Finn) Pint (0) + 2 0(F,Fp)ee (2) 


J=jo k 


[link] shows an example of a biorthogonal wavelet basis built by Cohen, Daubechies and Feauveau in [link], 
(called CDF-wavelets hereafter). 


phia4 phi.34 


psi34 psi.3¢ 


Primal and dual scaling and wavelet functions for the 
(3,1)-Cohen-Daubechies-Feauveau (CDF) 
biorthogonal basis. The primal wavelet function 


has one vanishing moment while the dual wavelet w 
has three vanishing moments. 


Fast wavelet transform 


One-dimensional wavelet transform 


Suppose we are given as signal the projection of a function onto the space 
Viq1: 
Equation: 


PFif= Yo spieGpie(®),  s441,k = (fF, bite) : 
k 


Using the dual refinement equations, we have: 
Equation: 


Sik = (fF, dix) = (f, So hidj+s2641 


= S Ai_2k$541,1 » 
k 


where the coefficients s ;, are called scaling coefficients, since they are 
related to scaling functions. Similarly, the wavelet or detail coefficientsd i, 
are obtained as 

Equation: 


dik = (f, dix) = S > G1-2n8 541, ; 
k 


The coefficients sj, and dj, are obtained from s;,17 by “moving average’ 
schemes, using the filter coefficients {hi} and {97} as ‘weights’, with the 


exception that these moving averages are sampled only at the even integers, 
i.e. a downsampling is performed. Such transform allows, once we have 


computed s 7 ,= ( f, b 1k) for a fine level J € N, to compute sj; and d jx 


for all coarser levels 7 < J without evaluating the integrals. 


Suppose now we are given the values of f at n = 27 equispaced design 
points. The scaling functions 7,4, k = 0,..., 2/7 — 1, are compactly 
supported and localized around 27k. Hence the coefficients ( hae b 1k) are 


weighted and scaled average of f on a neighborhood of 2-7 which 
becomes smaller as J tends to infinity. Consequently, it makes sense to 


replace the integral ( a4 db 14) by the (scaled) value of f at 2-7k. More 
complicate quadrature formulae have been developed in [link], [link], 
[link]. 

With 47 = {sack =0,...,2— 1} andd;:= {dajik=0,,...,2/ —1}, 
the forward (or analyzing) wavelet transform given by [link]-[link] can be 


rewritten as 
Equation: 


~ * ~* 
8; = Hjsj41 and dj =Gjsj41, 
~ * ~ 
where fH, denotes the Hermitian conjugate of ;. 
The inverse (or synthesis) transform is found by using the primal 


refinement equations and the fact that V;,; = V; @ Wj. 
Equation: 


Piif = > S5HAIP5+11 = »y Si,bP jk + S° dj,kW3,k 
k k 


> Sik » hi j+t.2K+1 + > dj k > GP j-41,2k-+1 
k l k I 
_ ss Pj+1,l (> hi—2k8 j,k + > a 
l k k 


| 


from which it follows that 
Equation: 


S411 = y hy_2K8 jk + S Gi—-2kd jx, - 
k k 


In matrix form, we have 
Equation: 


$741 = Hi;s;, ae G,d, * 


In the finite and classical setting, the matrices H,, G i ia ; and G ; are of 
size 27+! x 29. Moreover, if the basis functions are compactly supported, 


the four filters (hy, gi, hi, gi) have only a finite number of nonzero 
elements, and hence all these matrices are banded. 


Example: Haar wavelet transform 

~ * * 
In case of the orthogonal Haar transform, 4, = H, and is of the form 
Equation: 


ho hi 


ho hy 


since only hg and h, are different from zero : hy = hy = 1/ /2. The high- 


pass filter {g;} is such that g9 = —1/V2 and g; = 1/2. The forward 
transform [link]-[link] reduces to 
Equation: 


Sik = aD Sj4 e+ + 25 }+1,2k 
j gel; j 
J2 ie 
d 1 
ik — = $541,2k+1 — Ta $741,2k » 
J2 Fi 
and the reconstruction is given by 
Equation: 
1 1 d 
8; = — 
j+1,2k — 85k — — Us 
v2 v2 
1 1 
Sj+12k+1 = 85,4 + —=dyp 


V2 


Two-dimensional wavelet transform 


The wavelet transform has been successfully applied to compress images, 
which are modelled as functions defined on a regular two-dimensional grid. 


Sy Sy-1 


h 
Dy_4 


d 
DS 1 


Two-dimensional wavelet transform: first the filters are applied on the 
column of the matrix S'7, this produces two matrices. The filters are 
applied a second time on the columns of these two matrices, resulting 
in four elements: a matrix of scaling coefficients, and three detail 
matrices. 


The easiest way to build a two-dimensional MRA is probably to use tensor 
products of spaces, see [link], [link]. In terms of wavelet transforms, this 
leads to applying two times a one-dimensional transform: first on the ‘row' 
of the signal matrix S.z, and second on the “columns' of the resulting two 
matrices, see [link]. In this figure, we see that, at each level of the 
decomposition, three types of detail coefficients are produced: Di D; and 


Ds These superscripts recall that, in an image, horizontal edges will lead 
to large values of Dd. vertical edges will show up in D; and Dp will be 
sensitive to diagonal lines. 


However, such a transform is not able to compress efficiently an image that 
contains curves. More complex bidimensional bases are now proposed in 
the literature to better model discontinuities along curves, see for example 
[ink], [link], [link], [Link]. 


Non-decimated wavelet transform 


Suppose we have some signal {y;} observed at some equispaced design 
points : y; = f (i/n),i = 1,...,n withn = 27, J ¢N. The transform 
presented in the previous section is sometimes called ‘decimated’ because, 
for each scale 7, the coefficients dj, give only some information about the 
signal near the positions x = 2-/k, and not near all the existing design 
points 2-"k = k/n. 


For this reason, the decimated wavelet transform lacks the property of 
translation invariance: given to € R, the wavelet decomposition of f(. ) 
and of f(. —to) are in general completely different. This may lead to some 
unwanted pseudo-Gibbs oscillations near a discontinuity which is not 
localized at a dyadic point [link]. 


One remedy to this drawback consists in using a non-decimated wavelet 
transform (NDWT), also called translation-invariant (TT) [link] or 
stationary [link]. The idea behind the NDWT is to a perform a discrete 
wavelet transform, not only of the original sequence {y;};"_,, but of all the 
possible shifted sequences (Spy); = Y(t+h) mod n- In terms of wavelet 
functions, this transform corresponds to a set of functions 

Equation: 


Din (2) = (2) (@ —2°-%k)), J =Jo,--, J —1, F=0,...,27 -1. 


Ata given scale 7, the NDWT coefficients are thus present at all the 
locations k/n for k = 1,...,n and give information about the signal at each 
observed design point. In other words, the non-decimated transform fills in 
the gap introduced in the decimated transform, see [link]. 
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Schema illustrating the translation-invariant version of the Haar 
transform. The points marked by e are the one computed for the 
decimated Haar transform. At level J, one circulant shift is performed: 
the first observation is put at the end of the observed signal, and a 
second decimated transform is performed on the shifted data (yielding 
the points marked by o at level J — 1). This process is iterated at the 
coarser levels, producing detail coefficients at all the points. 


Since we have J — 79 scales and at each scales n detail coefficients, the 
NDWT gives an overdetermined representation of the original signal 
{y;};_, and the wavelet coefficients {dj,, 7 = 0,...,.J —1,k =1,...,n} 
are related to many different bases. Therefore the inverse operator will not 
be unique. A particular inverse, the average basis, corresponds to 
systematically average out the inverse wavelet transform obtained from 
each decimated wavelet transform that constitutes the translation-invariant 
transform. This makes the reconstruction robust with respect to a bad 
choice of a particular basis. Moreover, this average basis provides a 
smoother reconstruction than the original, decimated, transform [link], 
[link]. 


It allows for a (nearly) exact reconstruction of piecewise linear functions, 
instead of piecewise constant functions for the decimated Haar transform 
[link]. 


Approximation of Functions 


We first give a definition of the order of a multiresolution analysis. 


Note:(Order of a MRA in the classical setting) A multiresolution analysis is said to be 
of order NN if the primal scaling function ¢ reproduces polynomials up to degree 
N—1l,ie,For 0 <p < NV, ac, € such that «2? = S -cpe(a — fk) . 


The associated dual wavelet mn has then NV vanishing moments. In the classical setting, 
it is proved that the order of a MRA and the regularity of the scaling function are 
linked: the larger N, the higher the regularity of ¢. Symmetrically to [link], the order of 
the dual MRA is NN if db can reproduce polynomials up to degree N — 1. Figure 2 from 
Multiresolution analysis and wavelets shows an example of a biorthogonal basis where 


N = 3and N = 1. It illustrates the link between a high number of vanishing moments 
of the dual wavelet w and the regularity of the primal scaling function ¢. 


The main objective when decomposing a function in a wavelet series is to create a 
sparse representation of the function, that is, to obtain a decomposition where only a 
few number of detail coefficients are 

lar >!,whi < themaj or ityofthecoef ficientsareclose — zer © Bylarge', we 
mean that the absolute value of the detail coefficient is large. 


Near a singularity, large detail coefficients at different levels will be needed to recover 
the discontinuity. However, between points of singularity, we can hope to have small 


detail coefficients, in particular if the analyzing wavelets wp jk have a large number N of 
vanishing moments. Indeed, suppose the function f to be decomposed is analytic on the 


interval I without discontinuity. Since (a ; Bn) =O 10r pe = 0, sx, N — 1, weare 
sure that the first NV terms of a Taylor expansion of f will not give a contribution to the 
wavelet coefficient ( tf; Bn) provided that the support of Dik does not contain any 
singularities of the function f. 

This sparse representation explains why classical wavelets provide smoothness 
characterization of function spaces like the Hélder and Sobolev spaces [link], but also 


of more general Besov spaces, which may contain functions of inhomogeneous 
regularity [link], [link], [link], [link], [link]. 


We illustrate this characterization property with the case of 8—H6lder functions. 


Definition 2 
The class A® (L) of Hélder continuous functions is defined as follows: 


LifB <1, A9(L) = {f:1F(@) - Fl < De - ol? 

pean 
B>1,AP(L) = fF: [FD (a) — FD (y)| < Lie — yl? 5 |F2 (@)| <M}, 
where | (| is the largest integer less than 6 and 8’ = 6B — |B]. 


The global Hoélder regularity of a function can be characterized as follows [link], [link]. 


Note: Let f € A’ (L), and suppose that the (orthogonal) wavelet function w has r 
continuous derivatives and r vanishing moments with r > @. Then 


Equation: 


MPs djn)| < C2-HOHP) | 


A similar characterization exists for continuous and Sobolev functions [link], [link]. 


In the orthogonal setting, the wavelet ~ must be regular and have a high number of 
vanishing moments. On the contrary, in the biorthogonal expansion equation 5 from 


Multiresolution analysis and wavelets, it is mostly of interest to have a dual wavelet w 
with a high number of vanishing moments, and hence a regular primal scaling and 
wavelet functions. On the primal side, it is sufficient to have only one vanishing 
moment for wavelet denoising, and consequently wp may not be very regular. In this 


case, the wavelet coefficient ( f, Bn) with the less regular wavelet Wyk can be used to 
characterize f € A® (L) with0 < 6 < .N, evenif 8 > N = 1: witha biorthogonal 


basis, regular functions can be characterized by their inner products with much less 
regular functions. 


Nonparametric regression with wavelets 


In this section, we consider only real-valued wavelet functions that form an 


orthogonal basis, hence ¢ = b and = = w. We saw in Orthogonal Bases from 
Multiresolution analysis and wavelets how a given function belonging to 

Lz (IR) could be represented as a wavelet series. Here, we explain how to use 
a wavelet basis to construct a nonparametric estimator for the regression 
function m in the model 

Equation: 


Y; = m(az;) +6, i=1,...,n, n=2", JEN, 


where x; = + are equispaced design points and the errors are i.i.d. Gaussian, 
e~ N (0, 2), 


A wavelet estimator can be linear or nonlinear. The linear wavelet estimator 
proceeds by projecting the data onto a coarse level space. This estimator is of 
a kernel-type, see "Linear smoothing with wavelets". Another possibility for 
estimating m is to detect which detail coefficients convey the important 
information about the function m and to put equal to zero all the other 
coefficients. This yields a nonlinear wavelet estimator as described in 
"Nonlinear smoothing with wavelets". 


Linear smoothing with wavelets 


Suppose we are given data (x;, Y;);"_, coming from the model [link] and an 
orthogonal wavelet basis generated by {¢, =}. The linear wavelet estimator 
proceeds by choosing a cutting level 7; and represents an estimation of the 
projection of m onto the space V;j,: 

Equation: 


230-1 ji—123-1 


y= So 8inePjon (@) + >- ae dain (x) = S— 8j,nja (2), 
k=0 k 


J=Jo k=0 


with jo the coarsest level in the decomposition, and where the so-called 
empirical coefficients are computed as 
Equation: 


" <— ~ ins 
rn SCY; bj (@;) and dj, = - SCY; Wx (i) 
i=1 i=l 


The cutting level 7; plays the role of a smoothing parameter: a small value of 
41 means that many detail coefficients are left out, and this may lead to 
oversmoothing. On the other hand, if 7; is too large, too many coefficients 
will be kept, and some artificial bumps will probably remain in the estimation 
of m(z). 


To see that the estimator [link] is of a kernel-type, consider first the projection 


of m onto V;,: 
» ( / m (Yy) j,k widy) j,k (2) 


Equation: 
Py, m (zx) 
= [®: (x, y)m (y)dy , 


where the (convolution) kernel K;, (x, y) is given by 
Equation: 


=D, j,k ( YW) Pjr,k ( x) : 


Hardle et al. [link] studied the approximation properties of this projection 
operator. In order to estimate [link], Antoniadis et al. [link] proposed to take: 
Equation: 


—S_ il i/n 
Py,m(a) = Dv J Kn (muddy 
i=1 a n 


or 


i/n 


| 


Pji,k (y a) Pj,,k (x) . 


1)/ 


Approximating the last integral by +6,  (x;), we find back the estimator 
im (zx) in [linkl. 


By orthogonality of the wavelet transform and Parseval's equality, the L2— 
risk (or integrated mean square error IMSE) of a linear wavelet estimator is 
equal to the /2—risk of its wavelet coefficients: 


Equation: 
; ji-l . 9 
IMSE = E||m—mll7, = DBs 3x — Spal + D> Edu — 3,| 
=5o kk 
+ a = $1 + 52+ 83, 


JAN 


where 
Equation: 


Si := (mM, jz) and di = (m, jx) 


are called ‘theoretical’ coefficients in the regression context. The term 

S; + S¢ in [link] constitutes the stochastic bias whereas S3 is the 
deterministic bias. The optimal cutting level is such that these two bias are of 
the same order. If m is 6—H6lder continuous, it is easy to see that the optimal 
cutting level is j; (2) = O (n1/"+9)). The resulting optimal IMSE is of 


_ 28 | . . : 
ordern 7, In practice, cross-validation methods are often used to 
determine the optimal level 7; [link], [link]. 


Nonlinear smoothing with wavelets 


Hard-, soft-thresholding and wavelet estimator 


Given the regression model [link], we can decompose the empirical detail 
coefficient dip in [link] as 
Equation: 


n 
ma 


i 1 
din = = S- m, (a)P se (vi) + a S- E,W 5k (Li) 
= 


i=1 
= dik + Pik 


If the function m(z) allows for a sparse wavelet representation, only a few 
number of detail coefficients dj, contribute to the signal and are non- 


negligible. However, every empirical coefficient dj, has a non-zero 
contribution coming from the noise part p jx. 


Note: Note the link between the coefficients dz, in [link] and the theoretical 
coefficients diy, in [link]: 


Equation: 


dy, = —S~m(ai)bjn (2) 


i aoe ina: o(=) oe o(=) | 


In words, dj, constitutes a first order approximation (using the trapezoidal 
rule) of the integral d;,,. For the scaling coefficients s;,, it can be proved 


[link] that the order of accuracy of the trapezoidal rule is equal to N — 1, 
where WN is the order of the MRA associated to the scaling function. 


Suppose the noise level is not too high, so that the signal can be distinguished 
from the noise. Then, from the sparsity property of the wavelet, only the 
largest detail coefficients should be included in the wavelet estimator. Hence, 
when estimating an unknown function, it makes sense to include only those 
coefficients that are larger than some specified threshold value t: 

Equation: 


NH (dx, t) = dnl fig, os} 


This “keep-or-kill' operation is called hard thresholding, see [link |(a). 


Now, since each empirical coefficient consists of both a signal part and a 
noise part, it may be desirable to shrink even the coefficients that are larger 
than the threshold: 

Equation: 


diy = ns (dx,t) = sign (dj) (diet). . 


Since the function 7g is continuous in its first argument, this procedure is 
called soft thresholding. More complex thresholding schemes have been 
proposed in the literature [link], [link], [link]. They often appear as a 
compromise between soft and hard thresholding, see [link](b) for an example. 


In (a) the hard thresholding is represented in plain line: a coefficient dix 
with an absolute value below ¢ is put equal to zero. The soft thresholding 
is given in dashed line: there coefficients with absolute value above the 
threshold ¢ are shrunk of an amount equal to ¢. In (b), a more complex 
thresholding procedure, the SCAD threshold devised in Antoniadis and 
Fan [link] is represented. 


For a given threshold value ¢ and a thresholding scheme 77), the nonlinear 


wavelet estimator is given by 
Equation: 


fi (2) = So Bok bik (2) + Dome (Gast) dix (2) 
k j,k 


where jo denotes the primary resolution level. It indicates the level above 
which the detail coefficients are being manipulated. 


Let now d; = { din: k=0,...,27 — 1} denote the vector of empirical detail 


coefficients at level 7 and similarly define §,;. In practice a nonlinear wavelet 
estimator is obtained in three steps. 


1. Apply the analyzing (forward) wavelet transform on the observations 
{Y;};"_,, yielding s,, and d;, for 7 = jo,..., J — 1. 
2. Manipulate the detail coefficients above the level 7, e.g. by soft- 
thresholding them. 
3. Invert the wavelet transform and produce an estimation of m at the 
. a ~ n 
design points: {77 (x;)}._,. 
If necessary, a continuous estimator ™m can then be constructed by an 
appropriate interpolation of the estimated m (z;) values [link]. 


The choice of the primary resolution level in nonlinear wavelet estimation has 
the same importance as the choice of a particular kernel in local polynomial 
estimation, i.e., it is not the most important factor. It is common practice to 
take 79 = 2 or jo = 3, although a cross-validation determination is of course 
possible [link]. 


The selection of a threshold value is much more crucial. If it is chosen too 
large, the thresholding operation will kill too many coefficients. Too few 
coefficients will then be included in the reconstruction, resulting in an 
oversmoothed estimator. Conversely, a small threshold value will allow many 
coefficients to be included in the reconstruction, giving a rough, or 
undersmoothed estimator. A proper choice of the threshold involves thus a 
careful balance between smoothness and closeness of fit. 


In case of an orthogonal transform and i.i.d. white noise, the same threshold 
can be applied to all detail coefficients, since the errors in the wavelet domain 
are still i.i.d. white noise. However, if the errors are stationary but correlated, 
or if the transform is biorthogonal, a level-dependent threshold is necessary to 
obtain optimal results [link], [link]. Finally, in the irregular setting, a level and 
location dependent threshold must be utilized. 


Many efforts have been devoted to propose methods for selecting the 
threshold. We now review some of the procedures encountered in the 
literature. 


Choice of the threshold 


Universal threshold 


The most simple method to find a threshold whose value is supported by some 
statistical arguments, is probably to use the so-called ‘universal threshold’ 
[link], [link] 

Equation: 


tuniv = Og 2 log n , 


where the only quantity to be estimated is a. which constitutes the variance 
of the empirical wavelet coefficients. In case of an orthogonal transform, 


C4= 04) Ji 


In a wavelet transform, the detail coefficients at fine scales are, with a small 
fraction of exception, essentially pure noise. This is the reason why Donoho 
and Johnstone proposed in [link] to estimate og in a robust way using the 


median absolute deviation from the median (MAD) of d. Jui! 
Equation: 


median(|d jJ—1 — Median (4 7-1) ) 
0.6745 


Oqd = 


If the universal threshold is used in conjunction with soft thresholding, the 
resulting estimator possesses a noise-free property: with a high probability, an 
appropriate interpolation of {im (it produces an estimator which is at least 


as smooth as the function m, see Theorem 1.1 in [link]. Hence the 
reconstruction is of good visual quality, so that Donoho and Johnstone called 
the procedure ~VisuShrink' [link]. Although simple, this estimator enjoys a 
near-minimax adaptivity property, see ''Adaptivity of wavelet estimator". 
However, this near-optimality is an asymptotic one: for small sample size 
tuniv May be too large, leading to a poor mean square error. 


Oracle inequality 


Consider the soft-thresholded detail coefficients d’. Another approach to find 
an optimal threshold is to look at the J2—risk 
Equation: 


sh a ce dyn) = Bla dl, 


(3,k) 


and to relate this risk with the one of an ideal risk Aigeq). The ideal risk is the 
risk obtained if an oracle tells us exactly which coefficients to keep or to kill. 


In [link], Donoho and Johnstone showed that, when using the universal 
threshold, the following oracle inequality prevails 
Equation: 


o2 
& (@, d) < (2 log n+ 1) (<: + Fae : 
n 


However, this inequality is not optimal. Donoho and Johnstone looked for the 


optimal threshold t* (n) which leads to the smallest possible constant A. in 
place of 2 log n + 1. Such a threshold does not exist in closed form, but can 
be approximated numerically. For small sample size, it is sensibly smaller 
than the universal threshold. 


SureShrink procedure 


Given the expression [link] for the /2-risk, it is natural to look for a threshold 
that minimizes an estimation of this risk. 


By minimizing Stein's unbiased estimate of the risk [link] and using a soft 
thresholding scheme, the resulting estimator, called “SureShrink’, is adaptive 
over a wide range of function spaces including Hélder, Sobolev, and Besov 
spaces, see "Adaptivity of wavelet estimator". That is, without any a priori 
knowledge on the type or amount of regularity of the function, the SURE 


procedure nevertheless achieves the optimal rate of convergence that one 
could attain by knowing the regularity of the function. 


Other thresholding procedures 


We mention some of the other thresholding or shrinkage procedures proposed 
in the literature. 


Instead of considering each coefficient individually, Cai et al. [link], [link] 
consider blocks of empirical wavelet coefficients in order to make 
simultaneous shrinkage decisions about all coefficients within a block. 


Another fruitful idea is to use the Bayesian framework. There a prior 
distribution is imposed on the wavelet coefficients dj. This prior model is 
designed to capture the sparseness of the wavelet expansion. Next, the 
function is estimated by applying some Bayes rules on the resulting posterior 
distribution of the wavelet coefficients, see for example [link], [link], [link], 
[link]. 


Antoniadis and Fan [link] treat the problem of selecting the wavelet 
coefficients as a penalized least squares issue. Let W be the matrix of an 
orthogonal wavelet transform and Y := {Y;}/"_,. The detail coefficients 
d := {d;,} which minimize 

Equation: 


WY — dl, + > pa (Ide!) 
7k 


are used to estimate the true wavelet coefficients. In equation [link], p (-) is a 
penalty function which depends on the regularization parameter A. The 
authors provide a general framework, where different penalty functions 
correspond to different type of thresholding procedures (like, e.g., the soft- 
and hard- thresholding) and obtain oracle inequalities for a large class of 
penalty functions. 


Other methods include threshold selection by hypothesis testing [link], cross- 
validation [link], or generalized cross-validation [link], [link], which is used 
to estimated the /2-risk of the empirical detail coefficients. 


Linear versus nonlinear wavelet estimator 


In order to differenciate the behaviours of a linear and of a nonlinear wavelet 
estimator, we consider the Sobolev class Wj (C’) defined as 


Equation: 
We (C)=4 Ff: (FI E+ = 
#(C)= 4 FMR +l 


Fol so}, 


and that we denote V in short. Assume we know that m, the function to be 
estimated, belongs to V. In the next section, we will release this assumption. 
The L,,—risk of an arbitrary estimator T, based on the sample data is defined 
as E||T,, — ml|f, 1 < p < oo, whereas the L,—minimax risk is given by 
Equation: 


R,(V,p) =infsup E||T, — ml|>, 
In meV 


where the infimum is taken over all measurable estimators T;, of m. 
Similarly, we define the linear L,, —minimax risk as 
Equation: 


RY" (V,p) =infsup El|T;" — ml; 
a meV 


where the infimum is now taken over all linear estimators T™. Obviously, 
R™ (V,p) > Rn (V,p). We first state some definitions. 


Note: The sequences {a,,} and {b,,} are said to be asymptotically 
equivalent and are noted a, ~ b, if the ratio a,,/b, is bounded away from 


zero and co asn > co. 


Note:The sequence a,, is called optimal rate of convergence , (or minimax 


rate of convergence) on the class V for the L,,—risk if an ~ Rn(V, p) ue 
We say that an estimator m,, of m attains the optimal rate of convergence if 
suPmev E||mn — ml, ~ Rn (V,p). 


In order to fix the idea, we consider only the L2—risk in the remaining part of 
this section, thus p := 2. 


In [link], [link], the authors found that the optimal rate of convergence 
attainable by an estimator when the underlying function belongs to the 
Sobolev class W7 is an = n2-1, hence R, (V,2) = niet , We saw in 
"Linear smoothing with wavelets" that linear wavelet estimators attain the 
optimal rate for s—Hd6lder function in case of the D2—risk (also called 
‘IMSE'). For a Sobolev class W.’, the same result holds provided that q > 2. 


More precisely, we have the two following situations. 


1. If g > 2, we are in the so-called homogeneous zone. In this zone of 
spatial homogeneity, linear estimators can attain the optimal rate of 
convergence n~8/(28+1), 

2. If g < 2, we are in the non-homogeneous zone, where linear estimators 
do not attain the optimal rate of convergence. Instead, we have: 
Equation: 


Ri™ (V,2)/Rn(V,2) 3 00, as n > 00. 


The second result is due to the spatial variability of functions in Sobolev 
spaces with small index q. Linear estimators are based on the idea of spatial 
homogeneity of the function and hence do perform poorly in the presence of 
non-homogeneous functions. In contrast, even if q < 2, the SureShrink 
estimator attains the minimax rate [link]. The same type of results holds for 
more general Besov spaces, see for example [link], Chapter 10. 


Adaptivity of wavelet estimator 


We just saw that a nonlinear wavelet estimator is able to estimate in an 
optimal way functions of inhomogeneous regularity. However, it may not be 
sufficient to know that for m belonging to a given space, the estimator 
performs well. Indeed, in general we do not know which space the function 
belongs to. Hence it is of great interest to consider a scale of function classes 
and to look for an estimator that attains simultaneously the best rates of 
convergence across the whole scale. For example, the L,—Sobolev scale is a 
set of Sobolev function classes W,’ (C’) indexed by the parameters s and C, 
see [link] for the definition of a Sobolev class. We now formalize the notion 
of an adaptive estimator. 


Let A be a given set and let {.A., a € A} be the scale of functional classes 
F,, indexed by a € A. Denote R,, (a, p) the minimax risk over Fq for the 
L,—loss: 

Equation: 


R, (a,p) =inf sup Ellin — mll>. 
Mn MEF», 


* 
Note: The estimator m,, is called rate adaptive for the L,,—loss and the 
scale of classes Y,, a € A if for any a € A there exists c, > 0 such that 


Equation: 


sup E||m,, — mil, < Ca Rn (a,p) Vn > 1. 


MEF», 


The estimator m,, is called adaptive up to a logarithmic factor for the L,,— 
loss and the scale of classes Ay,a € A if for any a € A there exist cg > 0 
and Y = Ya > 0 such that 


Equation: 


sup Elm, — mil, < cqg(log n)"R, (a, p) Vn > 1. 


MCF » 


Thus, adaptive estimators have an optimal rate of convergence and behave as 
if they know in advance in which class the function to be estimated lies. 


The VisuShrink procedure is adaptive up to a logarithmic factor for the L2— 
loss over every Besov, Hélder and Sobolev class that is contained in C’'/0, 1], 
see Theorem 1.2 in [link]. The SureShrink estimator does better: it is adaptive 
for the L2—loss, for a large scale of Besov, Hélder and Sobolev classes, see 
Theorem 1 in [link]. 


Conclusion 


In this chapter, we saw the basic properties of standard wavelet theory and 
explained how these are related to the construction of wavelet regression 
estimators. 


